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ABSTRACT 

Bright circumstellar nebulae around massive stars are potentially useful to derive time-dependent 
mass-loss rates and hence constrain the evolution of the central stars. A key case in this context is the 
relatively young ejection-type nebula Ml-67 around the runaway Population I Wolf-Rayet star WR 124 
(=209 BAC), which exhibits a WN 8 spectrum. With WST-WFPC2 we have obtained a deep, Ha image 
of Ml-67. This image shows a wealth of complex detail that was briefly presented previously by Grosdi- 
dier et al. With the interferometer of the Universite Laval (Quebec, Canada), we have obtained comple- 
mentary Fabry-Perot Ha data using Canada-France-Hawaii Telescope (CFHT) MOS/SIS. From these 
data Ml-67 appears more or less as a spherical (or elliptical, with the major axis along the line of sight), 
thick, shell seen almost exactly along its direction of rapid spatial motion away from the observer in the 
ISM. However, a simple thick shell by itself would not explain the observed multiple radial velocities 
along the line of sight. This velocity dispersion leads one to consider Ml-67 as a thick accelerating shell. 

Given the extreme perturbations of the velocity field in Ml-67, it is virtually impossible to measure any 
systematic impact of the present WR (or previous LBV) wind on the nebular structure. The irregular 
nature of the velocity field is likely due to either large variations in the density distribution of the 
ambient ISM or large variations in the central star mass-loss history. In addition, either from the density 
field or the velocity field, we find no clear evidence for a bipolar outflow, as was claimed in other studies. 

On the deep Ha image we have performed continuous wavelet transforms to isolate stochastic structures 
of different characteristic size and look for scaling laws. Small-scale wavelet coefficients show that the 
density field of Ml-67 is remarkably structured in chaotically (or possibly radially) oriented filaments 
everywhere in the nebula. We draw attention to a short, marginally inertial range at the smallest scales 
(6.7-15.0 x 10 -3 pc), which can be attributed to turbulence in the nebula, and a strong scale break at 
larger scales. Examination of the structure functions for different orders shows that the turbulent regime 
may be intermittent. Using our Fabry-Perot interferograms, we also present an investigation of the sta- 
tistical properties of fluctuating gas motions using structure functions traced by Ha emission-line cen- 
troid velocities. We find that there is a clear correlation at scales 0.02-0.22 pc between the mean 
quadratic differences of radial velocities and distance over the surface of the nebula. This implies that the 
velocity field shows an inertial range likely related to turbulence, though not coincident with the small 
inertial range detected from the density field. The first- and second-order moments of the velocity 
increments are found to scale as < | Av(r) \ ) ~ r 0,5 and < | Av(r) ] 2 3 > ~ r 0,9 . The former scaling law strongly 
suggests that supersonic, compressible turbulence is at play in the nebula; on the other hand, the latter 
scaling law agrees very well with Larson-type laws for velocity turbulence. Examination of the structure 
functions for different orders shows that the turbulent regime is slightly intermittent and highly multi- 
fractal with universal multifractal indexes a « 1.90-1.92 and C, « 0.04 + 0.01. 

Subject headings: ISM: bubbles — ISM: individual (Ml-67) — ISM: kinematics and dynamics — 
stars: individual (WR 124) — stars: Wolf-Rayet — turbulence 
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1. INTRODUCTION 

1.1. Nebulae Ejected by Hot Stars 
Nebulae around central stars of planetary nebulae 
(CSPN) or Population I Wolf-Rayet (WR) stars have 
proven to be a powerful diagnostic tool to understand how 
mass is lost in the post main-sequence phases of stellar 
evolution. CSPN and some («j; Marston 1995, 1996, 1999) 
massive WR stars are surrounded by circumstellar nebulae 
(CN) made up of stellar material ejected during successive 
stellar wind phases. On the whole, these nebulae often 
appear to have similar physical properties, which point to a 
common formation mechanism (Chu 1993) independent of 
the strong differences between both types of hot stars. 

The dynamics of planetary nebulae and CN surrounding 
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Population I WR stars are extremely sensitive to the stellar 
wind velocity and mass-loss history. In the context of the 
interacting stellar wind model, a spherically symmetric, fast, 
hot stellar wind (r cxp ~ 10 3 km s _1 ) is catching up and 
colliding with a preexisting slower (r cxp ~ 10 1 —1 0 2 km s _1 ), 
denser wind (CSPN: Kwok, Purton, & FitzGerald 1978; 
Chu et al. 1991; Balick 1994; Population I WR stars: 
Marston 1995, 1996, 1999; Garcia-Segura & Mac Low 
1995; Garcia-Segura, Mac Low, & Langer 1996, and refer- 
ences therein). The noncircular shapes of many CN suggest 
that the fast wind is sweeping up a slow wind that is flat- 
tened by rotation or binary effects (Garcia-Segura & Mac 
Low 1995; Mellema & Frank 1995). However, a flattened, 
fast wind blowing out a spherically symmetric slow wind 
can also lead to a noncircular CN (Frank, Ryu, & Davidson 
1998). 

Rayleigh-Taylor instabilities are expected to occur as a 
consequence of such interaction, as well as other hydrody- 
namical instabilities that may lead to turbulence in the 
nebula. On a qualitative basis, the interacting stellar wind 
model agrees with the observations, but it is still difficult to 
determine the precise initial, overall geometry of the collid- 
ing winds. In particular, the question arises as to the impact 
of precollision, apparently universal, clumping or possible 
turbulence of the hot-star winds (Robert 1992; Balick et al. 
1996; Lepine, Moffat, & Henriksen 1996; Acker, Grosdi- 
dier, & Durand 1997; Grosdidier et al. 1997; Eversberg, 
Lepine, & Moffat 1998; Lepine & Moffat 1999; Grosdidier, 
Acker, & Moffat 2000, 2001 ; and references therein) on the 
dynamics/morphology of CN. Moreover, despite prelimi- 
nary work done by Stone, Xu, & Mundy (1995), the inci- 
dence of a highly variable (as well in the spatial as in the 
time domain) mechanical luminosity [ oc M(r, 9, <f), t) 
x Vf lnw (r, 9, (j>, tj] still has to be investigated in detail. In 
particular, the turbulence of the ejected nebulae has to be 
studied in detail in order to test whether it is related to 
preexisting turbulence (or any other process of fragmenta- 
tion, hence variability) that already arose in the winds them- 
selves before they interacted. 

Recall that the spectral variability originating in massive 
WR stars and CSPN shows the WR phenomenon to be 
remarkably similar in both (Grosdidier et al. 2000, 2001). 
Despite the strong differences between both types of hot 
stars, this supports the understanding of the WR spectral 
phenomenon as being purely atmospheric. This, in line with 
the believed similarity of the processes that are at play in the 
formation of the related CN (Chu 1993), places one in a 
position to check also for universality of the shaping of CN 
surrounding any hot stars. 

1.2. Turbulence in Ejected Nebulae and H n Regions 

Like H n regions, gas motions in nebulae ejected by hot 
stars are characterized by very large Reynolds numbers. 
For example, the Reynolds numbers exhibited by typical H 
ii regions may be as large as 10 9 (Boily 1993), using 1 pc for 
the characteristic spatial scale, 10 km s _1 for the velocity 
scale, and considering the normal kinematic viscosity for 
ionized hydrogen gas at electronic temperatures of about 
10,000 K. Such values are well above the critical values 
(10 2 -10 4 ) determining the transition to turbulence. Like H ii 
regions, it is expected that the kinematics of nebulae ejected 
by hot stars would be best described in terms of turbulent 
flows. Turbulence is characterized by a high degree of non- 
linearity in the governing equations and complex coupling 


mechanisms, leading to a variety of physical/transport pro- 
cesses occurring over a great variety of scale lengths. 

Turbulence is inherently apparently stochastic, since the 
flow variables (e.g., velocity, pressure, density) fluctuate in 
an unpredictable manner about their mean values. 
However, turbulence also exhibits specific length scales and 
scaling laws, with rapidly increasing numbers of smaller 
features as a result of energy cascading down a ladder 
leading ultimately to dissipation (Fleck 1996). At small 
scales, turbulent energy returns to the ambient medium as 
thermal energy. The nearly dissipationless energy cascade 
makes the flow not strictly unpredictable if one attempts to 
study quantitatively its structure via a statistical approach. 
Indeed, the signature of turbulence is confirmed by the pres- 
ence of correlations in the flow variables between neighbor- 
ing points that can be detected by means of statistical 
functions (Scalo 1984). The nonlinearity of the equations of 
hydrodynamics causes high-amplitude structures to change 
their shape in a way that continuously broadens the wave 
spectrum to include smaller and smaller scales. A simple, 
unique wave train is expected to steepen into a shock front 
damping at the leading edge. However, an intricate pattern 
of motions likely drives a cascade of energy dissipation 
taking on a fractal structure (Henriksen 1994 and references 
therein). 

A statistical search for any degree of correlation between 
neighboring points could be carried out via two main 
approaches: autocorrelation functions or structure func- 
tions (these two statistical tools are more reliable than the 
dispersion-scale relation); Scalo 1984; Miville-Deschenes & 
Joncas 1995). The so-called structure function is defined as 
the average of the quadratic differences of velocities as a 
function of the spatial separation between the points where 
the velocity is measured (Miesch & Bally 1994). In this 
study we favor the structure function calculation because it 
is more robust than the autocorrelation method. The 
robustness of the structure function technique originates in 
one important fact : its application requires the use of veloc- 
ity differences rather than velocity “ coherence ” (in the case 
of the autocorrelation method) between points separated by 
a given distance. Indeed, the use of the autocorrelation func- 
tion assumes that the statistical properties of the flows are 
not a function of position, the velocity and density fields 
being considered stationary. When using structure func- 
tions, we assume a less stringent hypothesis: local station- 
arity. However, note that at large scales the flows always 
lose their stationarity (Scalo 1984). The aim of the present 
paper is to describe the turbulent behavior of nebulae 
ejected by hot stars, both from the density structure (via 
high-resolution imagery) and the kinematical structure (via 
high spatial resolution Fabry-Perot interferometry). The 
search for turbulence will be carried out through the calcu- 
lation of the moments of the density/velocity increments for 
different orders (i.e., we will not limit our investigation to 
the second-order structure function, as is generally the case 
in other studies). Examination of the structure function for 
different orders additionally provides a way to estimate the 
level of intermittency in the flows. 

Ml-67, surrounding the WR star WR 124 (WN 8), was 
chosen for its relatively large apparent angular diameter 
(which allows one to study any possible energy cascade 
within a large span of scale lengths) and its suspected weak 
interaction with the ambient interstellar medium. The latter 
condition ensures that (1) the ejecta would be less affected 
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Fig. 1. — Deep WFPC2/Ha log-scale image of Ml-67 from four combined exposures totalling nearly 3 hr. The field stars have been subtracted out 
(enlargement of Fig. 1 from Grosdidier et al. 1998). 


by collision and mixing with the preexisting ambient gas, 
and therefore (2) the dynamics and structure of the ejected 
nebula would best reflect the initial conditions of the ejected 
winds. This, in line with the apparently universal, ubiqui- 
tous clumping of hot stellar winds, places one in a position 
to test the possible impact of the wind fragmentation on the 
nebular fragmentation. Note that high resolution optical 
spectra of WR 124 recently obtained by A. Acker (2000, 
private communication) at the Observatoire de Haute 
Provence (France) and Marchenko et al. (1998) results 
indeed show the central star of Ml-67 to be variable. 

2. THE NEBULA Ml-67 AND ITS CENTRAL STAR WR 124 

Ground-based, narrowband images and discussion of the 
ejected nebulosity Ml-67 surrounding the cool nitrogen 


WN 8 WR star WR 124 ( = 209 BAC) are given in Esteban 
et al. (1991, 1993), Sirianni et al. (1998), Crowther et al. 
(1999), and references therein. Until 1991 the central star of 
Ml-67 was often mistaken for a CSPN. From a study of the 
Na i D 2 interstellar absorption line in the context of Galac- 
tic rotation, Crawford & Barlow (1991) showed that the 
distance to WR 124 is 4-5 kpc. For such a large distance, 
the central star must be a massive WR star rather than a 
CSPN, given its relatively high apparent brightness, 
V x 11.6 (Massey 1984), and the huge interstellar reddening 
to Ml-67 ( E b _ v x 0.9-1. 5; Solf & Carsenty 1982; Esteban 
et al. 1991; Crowther et al. 1995b). This is compatible with 
M v x — 6.0, which is typical for WN 8 stars. Photoioniza- 
tion modeling of Ml-67 has been made by Crowther et al. 
(1999), who were able to trace the central star H-Lyman 
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Fig. 2. — Enlargement of the faint outer northeastern region (P.A. « 10°-30°) of the nebula ( upper panel) where at least nine “reversed” bow shocks are 
visible ( lower panel). Such features are seen throughout the nebula although better detected at the periphery. A few residuals of subtracted field stars also 
appear. The chosen contrast emphasizes the dimmer features. 


continuum flux distribution, which is normally not directly 
detectable. The main outcome of this work was to establish 
the crucial importance of a line-blanketed energy distribu- 
tion for proper modeling of the WN 8 atmosphere. 

As an E-type (ejected material) nebula (see Chu, Treffers, 
& Kwitter 1983; Smith 1995 and references therein), Ml-67 
constitutes an object in its earliest phase of wind inter- 
action, when the ejecta are least affected by collision and 
mixing with the interstellar medium, and thus best reflect 
the initial conditions of the ejected winds. This is supported 
by detailed abundance analysis of the nebula: N/O « 2.95 
(Esteban et al. 1991), rather than N/O* 0.07 for the Galac- 
tic ISM (Shaver et al. 1983). Note that the E-type Galactic 
CN RCW 58 (central star: WR 40; spectral type: WN 8), is 
also obviously enriched and must contain stellar ejecta; 
however, this nebula already shows a circumstellar ring, 
which suggests a more pronounced interaction with its 
ambient medium. For the same reason, we do not consider 
the multiple ring (Marston 1995) CN surrounding the WN 
8 star WR 16: the two outer rings are likely E-type ring 
nebulae, but their dynamical ages are relatively long (some 
10 6 yr and 7 x 10 s yr). In addition, the inner ring is likely a 
W-type (wind-blown bubble; see Chu et al. 1983 and refer- 
ences therein) CN. 


Ground-based observations (e.g., Chu & Treffers 1981; 
Sirianni et al. 1998) exhibit a nearly spherical, possibly 
bipolar, patchy structure made up of many bright unre- 
solved knots and filaments. The basic properties of Ml-67 
are: diameter 110"-120" (i.e., 2.4— 2.6 pc at a distance of 4.5 
kpc; Grosdidier et al. 1998), mean expansion velocity 40-45 
km s -1 (Sirianni et al. 1998), implying a dynamical age of 
the order of a few 10 4 yr. 

In the framework of our search for the direct influence of 
a clumpy, turbulent-like stellar wind on the surrounding 
nebula it ejects and more generally on the interstellar 
medium (Moffat 1998), we have obtained Hubble Space 
Telescope (HST) Wide Field Planetary Camera Two 
(WFPC2) images of the nebula Ml-67. Figure 1 is an 
enlargement of the picture reproduced in Grosdidier et al. 
(1998). They reported no overall global shell structure to the 
nebula (see also Grosdidier et al. 1999) from the projected 
radial distribution of the azimuthally averaged Ha surface 
brightness. Rather, the radial profile centered on WR 124 
agrees with a decreasing, wind-density power-law distribu- 
tion with increasing distance from the central star. Grosdi- 
dier et al. (1998) understood the absence of a well-defined 
shell as evidence for an extremely small age of the nebula : 
the WR phase may have just turned on, making the inter- 
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action zone (i.e., the nebula) of the present fast wind with the 
previous, slower wind still undetectable (which is in line 
with the presence of El emission lines in the spectrum of WR 
124; Crowther et al. 1995a). 

However, these deep HST/WFPC2 images in Ha of the 
nebula also suggested that Ml-67 may likely be the imprint 
of a previous, luminous blue variable (LBV) slow wind 
ejected from the WR central star WR 124, which is sup- 
ported by other arguments (Sirianni et al. 1998, and refer- 
ences therein; Massey & Johnson 1998) as opposed to a 
classical red supergiant very slow wind (Smith 1996) blown 
and compressed by the present WR wind. Note that ubiqui- 
tous reversed bow-shock-like structures are detected 
throughout the nebula (the more prominent ones being 
detected at the outer edge of Ml-67, P.A. « 10°-20° and 
130°-160°; see Fig. 2). Generally, the whole periphery of the 
nebula appears made up of numerous small (l"-5") reversed 
bow-shock-like structures. These features can be inter- 
preted as the result of the thin-shell instability of the highly 
variable previous LBV stellar wind interacting with itself 
(e.g., Stone et al. 1995). In addition, some dense, persisting 
clumps have possibly been ejected directly from the LBV 
stellar surface (Grosdidier et al. 1998). 

3. OBSERVATIONS AND DATA REDUCTION 
3.1. HST Imagery 

Here, we expand on the description given in Grosdidier et 
al. (1998). Four WFPC2 images of Ml-67 with a total com- 
bined exposure of 10,000 s were taken in 1997 March with 
the HST through the narrowband F656N Ha filter (Biretta 
et al. 1996). They were combined (with the task GCOM- 
BINE in STSDAS/IRAF 9 ) to produce a single deep Ha 
image cleaned from bad/hot pixels and cosmic rays. Before 
the combination of the four F656N Ha images, we applied 
proper shifts to them in order to achieve alignments to 
within + 0.01-0.03 pixels. The shifts were determined using 
at least 15 (Planetary Camera; PC) and 25-30 (Wide Field 
Cameras 1, 2, and 3; WFs) bright stars and calculating the 
relative displacements between each exposure through the 
estimation of the star centroids. Note that rather than 
applying one shift to an image to match the other, we 
applied half-displacements to both images. This reduces the 
differences in artificial broadening of the stellar profiles 
between the two images. Finally, we achieved full widths at 
half maximum spanning 2.4-2.5 pixels ( « 0" 1 1) for the PC, 
and 1. 7-2.0 pixels (0"17-0"2) for the WFs, in the Ha band. 

The four narrowband F656N Ha images were straddled 
in wavelength by broadband V and R images : four WFPC2 
images of the same field taken through each of the broad- 
band F675W R and F555W V filters (Biretta et al. 1996) 
were separately combined in the same way to produce 
similar high signal-to-noise ratio images of “ continuum ” 
starlight relatively close to the wavelength of Ha. The V and 
R images enabled us to interpolate the stellar light to the 
wavelength corresponding to Ha. After proper flux scaling 
and spatial shifting (also within +0.01-0.03 pixels), the 
interpolated broadband image 10 was carefully subtracted 

9 IRAF is distributed by the National Optical Astronomy Observa- 
tories, operated by the Association of Universities for Research in 
Astronomy, Inc., under cooperative agreement with the National Science 
Foundation. 

10 The F-band component being almost negligible. 


from the Ha image in order to produce a deep continuum- 
subtracted Ha image with the field stars removed. Note that 
we achieved full widths at half maximum spanning 2.4-2.5 
pixels for the PC, and 1. 7-2.0 pixels for the WFs, both in the 
R and V bands. These values are similar to those obtained 
for the combined Ha image and permitted an optimum 
continuum starlight subtraction. 

The scaling in intensity between the interpolated broad- 
band image and the combined Ha image was made by mea- 
suring the flux of at least 20 bright stars in each chip, and 
minimizing the residuals in the continuum-subtracted Ha 
image. This procedure works quite well for Ha line-free 
stars, which are the majority. For WR 124 we applied a 
different flux scaling because this WN 8 star emits a copious 
amount of Ha + He n 26560 light. It is worthy of note that 
we have obtained all images within a contiguous span of 
time (six HST orbits). This reduced problems with potential 
variable stars in the subtraction of starlight from the orig- 
inal Ha image. Note that the images were originally pipeline 
processed before being released but required a full recalib- 
ration because of severe changes in the dark files. 

3.2. Fabry-Perot Interferometry 

Nebular kinematics can be studied very effectively by 
observations at optical wavelengths with the use of a scan- 
ning Fabry-Perot interferometer in combination with a 
two-dimensional detector. This technique allows one to 
probe the velocity field with a wider field of view and higher 
throughput than with a normal long-slit, high-dispersion 
spectrograph. 

With the servo-stabilized Fabry-Perot (FP) interferome- 
ter of the Universite Laval (Joncas & Roy 1984), we 
obtained complementary FP Ha data of Ml-67 using 
MOS/SIS 11 (Le Fevre et al. 1994) at the Canada-France- 
Hawaii Telescope (CFHT), in 1996 August 30. The instru- 
ment was mounted at the Cassegrain f/8 focus and covered 
a field of about 200" on a side. An Ha narrowband filter 

o o 

(<S2 = 8.8 A) centered at 6575 A was used as a premonochro- 
mator. The combination of tilting and air temperature blue- 
shifted the central wavelength of the filter to match the 
red-shifted Ha line originating in Ml-67 (mean peculiar 
velocity: about + 137 km s” 1 ; Sirianni et al. 1998). The FP 
interferometer has a finesse F = 30 and a free spectral range 
of 392 km s -1 at Ha, where the interference order is 
p = 765. Since both the interferometer and filter were tilted, 
we avoided ghost images (Georgelin 1970). We used the 
2048 x 2048 pixel STIS2 CCD (gain = 4 e~ ADU 1 , 
readout noise = 8 e~, pixel width = 21 g) with a spatial 
sampling of 0'.'30 per pixel. We employed the MOS/SIS 
tip/tilt image stabilizer in order to improve image quality 
and achieved effective seeing spanning the range 
FWHM « 0'.'6-0'.'7. The useful data cube was 1500 x 1500 
pixels in the spatial domain, and 66 channels (2.2 x F, in 
order to satisfy the Nyquist criterium) in the wavelength 
domain, giving a velocity sampling of 5.9 km s -1 per 
channel. The integration time was 200 s channel -1 . A neon 
calibration lamp was used to obtain rest-frame interfero- 
grams. Finally, the data processing was done using the 
IRAF and ADHOC 12 (see Amram et al. 1995, Blais- 
Ouellette et al. 1999) packages. 

1 1 http ://www.cfht.hawaii.edu/Instruments/Spectroscopy/SIS. 

12 http://www-obs.cnrs-mrs.fr/ADHOC/adhoc.html. 
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4. THE TURBULENT STATUS OF Ml-67 

4.1. HST Imagery: Wavelets for Structure Function 
Analysis 

4.1.1. Method 

We favor the use of wavelets to perform structure func- 
tion calculations because this method is much more robust 
than the classical study of the scaling behavior via the 
power spectrum estimation as a function of the modulus of 
the wavevector k. Indeed, the latter method is sensitive to 
biases in the data, e.g., smooth polynomial trends, that can 
easily be shown to drastically alfect the power spectrum 
scaling law. In addition, this method is appropriate when 
the fractal signal under consideration is not associated to a 
unique roughness exponent H (Schmittbuhl, Vilotte, & 
Roux 1995). 

Using two-dimensional wavelets, we have analyzed our 
deep continuum-subtracted Ha image of Ml-67 to isolate 
stochastic structures of different characteristic size. Hence 
we effectively looked for scaling laws, as has been clearly 
demonstrated in other astrophysical contexts (e.g., clump- 
ing in CO: Gill & Henriksen 1990; clumping in WR winds: 
Lepine 1994). In what follows, the deep Ha image will be 
considered as a two-dimensional function noted /. 

Following Muzy, Bacry, & Arneodo (1993), we have 
explicitly used wavelets to perform a structure function 
analysis. We have considered the field of “ increments ” of / 
over a projected spatial distance r at the pixel (x, y), 

A flrix, >>)] =/[(*, y) + R/2] -/[(*, y) - R/2] , (1) 

by the continuous wavelet transform of / at the scale r, 
T w (f)[r;(x, y)], 

A flrix, y)] « T v (f)[r;(x, y)] 

= r~ 2 x (/*'F r )[r;(x, y)] , (2) 

where R is a two-dimensional increment vector (on the 
plane of the sky) of length r, * denotes a convolution, and 
V F,. is the analyzing wavelet at the spatial scale r (Muzy et al. 
1993). In our study, the so-called mother-wavelet, V F, is the 
two-dimensional “Mexican hat” (formally the second 
spatial derivative of e~ ix2+y2)l2 : (2 — x 2 — y 2 )e~ <x2+y2)/2 ) 
and v P r (A, y) = v P(x/r, y/r) (Farge 1992). It is worth noting 
that the Mexican hat wavelet is orthogonal to polynomials 
of orders 0 and 1. Therefore, the wavelet transform is insen- 
sitive to linear trends, hence linear nonstationary com- 
ponents, in the signal. 

The statistical moments of | A f[r;(x, y)] | (the two-point 
correlation known as “ structure function of order p ”) were 
estimated through the statistical moments of 
I T,,( /j[r ;(x, y)] | by averaging over the location : 

<| A/(r)|"> * cj | Uf)[r-(x, y)] fdxdy , (3) 

where C is an arbitrary constant. Note that Muzy et al. 
(1993) applied this formalism to one-dimensional signals. 
Here we generalize their approach to two-dimensional 
signals. In § 4.1.2 we will perform tests on artificial fractal 
two-dimensional signals and show that this procedure is 
correct. Since this approach does not depend on the analyz- 
ing wavelet (Muzy et al. 1993), one is in a position to test for 
scaling laws of the form: < | A/(r) | p > ~ r C(p) . 

Note that ((0) = 0 by definition and ifp) is a smooth, 
differentiable, monotonically nondecreasing, concave func- 


tion of p (if the signal has absolute bounds), no matter how 
rough the data are (e.g., Marshak et al. 1994). Continuous 
signals satisfy < | A/(r) | ) ~ r, hence C{p) = p. Processes with 
tfp) oz p are called “ monoaffine ” or “ monofractal,” 
whereas processes with variable C(p)/p are called 
multiaffine ” or “ multifractal.” 

4.1.2. Tests on Artificial Two-dimensional Fractal Signals 

In order to verify the validity of equation (3), it is neces- 
sary to test our numerical code on artificial two- 
dimensional signals. Preferably this should be done on 
two-dimensional fractals in order to check the efficiency of 
our method in finding self-similarity, hence fractality, in the 
fields of increments. 

The basic method for creating artificial fractals is to con- 
sider the spectral density, that is, the mean square fluctua- 
tion at any particular frequency and how that varies with 
frequency. This technique, the power spectral method, is 
based upon the observation that many natural forms and 
signals have a l/v p frequency power spectrum, i.e., the spec- 
trum falls off as the inverse of some power of the frequency 
v, where the spectral component (I is related to the fractal 
dimension of the signal: ft = 2 H + 2, with the persistence 
parameter 0 < H < 1. The latter parameter quantifies the 
degree of “ roughness ” (e.g., H = 0 for a fully rough signal; 
H = 1 for a fully smooth signal) and is related to the fractal 
dimension, D : H = 3 — D (Moghaddam, Hintz, & Stewart 
1991; Cox & Wang 1993). Thus, D will take on values 
between 2 and 3. 

Figures 3, 4, and 5 show three constructions of artificial 
fractal signals along with a horizontal cut of the images. To 
build these signals we first generate a random 512 x 512 
white-noise signal. We then apply a Fourier transform to 
this image leading to a secondary signal, which essentially 
looks like another noise field. We scale the resulting spectra 
by the desired l/v f function and finally apply an inverse 
Fourier transform. Controlling how rough or smooth the 
surface is has to be done by varying the power relationship. 
From Figure 3 to Figure 5 we have increased the value of ft, 
hence decreasing the roughness. Note that these three 
images are all based on the same initial seed and therefore 
have the same general shape. 

Figure 6 (left panels) shows our numerical applications of 
equation (3) for p = 1 to the three artificial fractal signals. 
Note that choosing p = 1 in equation (3) and the r " 2 scaling 
in equation (2) relates any power law < | A/(r) | ) ~ r C(1) to 
the persistence parameter H : H = ((1). Inspection of Figure 
6 suggests that the conjecture given in equation (3) is very 
accurate: the derived slopes are always very close to the 
theoretical value of H. 

The right panels of Figure 6 show the hierarchy of expo- 
nents C(p) related to the signals shown in Figures 3, 4, and 5. 
Here the statistical moments of | A/' [r ;(x, y )] | have been 
estimated for each artificial signal, for p = 1-5, in steps of 
0.5. The search for scaling laws has been performed as 
before for any individual value of p. Note that our artificial 
signals are monofractal, therefore £(p) should be strictly 
proportional to p. This is seen very clearly for H = 0.05. 
However, for H = 0.475 and H = 0.9 we detect a clear 
departure from this theoretical property for p > 3 and 
p > 1, respectively. This is mainly due to the small ampli- 
tude fluctuations of the last two signals (see Figs. 4 and 5, 
lower panels) : for large values of p, application of equation 
(3) may lead to very small numbers whenever 
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Fig. 3. — Artificial fractal signal with H = 0.05 (upper panel), along with 
a horizontal cut in the middle of the image ( lower panel). 




Fig. 4. — Artificial fractal signal with H = 0.475 ( upper panel), along 
with a horizontal cut in the middle of the image ( lower panel). 


I T,,(/')[r ;(x, >’)] | < 1 , in combination with dominant 
extreme events, which are poorly sampled by definition. 

These numerical experiments demonstrate how our 
wavelet method is efficient in detecting any fractal structure. 
However, the reliability of the structure function analysis 
for the highest orders appears very sensitive to the ampli- 
tude of the fluctuations. 

4.1.3. Results on the HST /Ha Deep Image of Ml-67 

Figures 7, 8, and 9 show, respectively, the wavelet coeffi- 
cients of the deep, net Ha image (see Fig. 1) for three partic- 
ular spatial scales. At the largest scale, 11 "3 (Fig. 7), the 
wavelet coefficients reveal the large-scale structure of 
Ml-67, where one can easily recognize the bright arcs of 
Figure 1. The number of clumps then increases rapidly 
towards smaller scales (see Figs. 8 and 9). More strikingly, 
at the smallest scale close to the resolution limit, 0'.'3 (Fig. 9), 
the density field of Ml-67 appears remarkably structured in 
apparently chaotically oriented filaments everywhere in the 
nebula. However, many filaments (especially at the smallest 
scales) seem to be slightly radially oriented. Since the field 
stars have been subtracted off in Figure 1, these filaments 


are strictly related to the fluctuating density field of Ml-67. 
Note that neither at large scales (Fig. 7), nor at small scales 
(Fig. 9) is any bipolar outflow clearly detected. Using the Ha 
HST image, Figure 10 shows the projected angular dis- 
tribution of the surface brightness averaged over 10°-wide 
sectors centered on WR 124. From the median two possible 
lobes are detected for P.A. « 65°, and P.A. « 185°. Since 
the lobes are not antialigned on the plane of the sky, they 
are unlikely related to a bipolar outflow. Considering the 
mean, the situation is even worse: many subpeaks reveal the 
highly irregular and patchy structure of Ml-67. Using wider 
sectors (hence losing angular resolution), we arrive at the 
same conclusions. 

In Figure 11 the quantities < | A/(r)| p > are plotted down 
to 3 pixels (0'.'3) of WFPC2, in order to secure good sam- 
pling of the analyzing wavelet, for p = 1 and 2. (The contri- 
bution of the instrumental noise has been subtracted off by 
performing the same analysis on a rectangular region 
outside the nebula; the field stars do not contribute in this 
calculation since they have been removed.) Note the short, 
marginal inertial range (less than a decade) at the smallest 
scales from about 6.7 x 10“ I * 3 to about 1.5 x 10 -2 pc, which 
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D=2.10 



Fig. 5. — Artificial fractal signal with H = 0.90 (upper panel), along with 
a horizontal cut in the middle of the image ( lower panel). 


could be attributed to turbulence. More strikingly, we 
detect the presence of a strong scale-break at scales of 
<L> « 5"-10" (or 0.11-0.22 pc, adopting a distance of 4.5 
kpc for Ml -67), corresponding to the characteristic widths 
of the bright arcs seen in Figure 1. This suggests that the 
bright arcs may be related to a specific earlier ejection event 
or hydrodynamic instability. Looking carefully at Figure 11 
a second characteristic scale is found at <s> « l"-2" (or 
around 3.3 x 10” 2 pc), where the second-order structure 
function is still growing but at a significantly lower rate. We 
interpret these characteristic scales as the presence of essen- 
tially two types of structures in the signal: bright, small 
<s>-scale features, and relatively less intense, larger <L>- 
scale features. Indeed, as seen in the second-order structure 
function which emphasizes the most intense events, as the 
scale r increases <|A/(r)| 2 > first increases sharply because 
the increments from the bottom to the top of the <s>-scale 
features dominate and there are more and more of these. 
Then <|A/(r)| 2 > becomes saturated with respect to these 
features: for all scales above <s>, <|A/(r)| 2 > contributes a 
relatively fixed number of events to the spatial average. 
Further on we see the larger <L>-scale features at work in a 



Fig. 6 . — Left panels: Persistence parameter retrieval for H = 0.05, 
H = 0.475, and H = 0.9 (cf. Figs. 3, 4, and 5). Right panels: Structure 
function analysis of the signals shown in Figs. 3, 4, and 5. 



Fig. 7. — Wavelet coefficients of Ml-67 at the scale 1173. The orienta- 
tion of the image is the same as of Fig. 1. 
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Fig. 8. — Wavelet coefficients of Ml-67 at the scale 1"8. The orientation 
of the image is the same as of Fig. 1. 
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Fig. 10. — Projected angular distribution of the surface brightness aver- 
aged over 10°-wide sectors centered on WR 124. 


similar way. They allow < | A/(r) | 2 > to keep growing but at 
a slower rate because smaller increments are spawned by 
these features. For scales larger than <L>, the second-order 
structure function again becomes saturated. For scales close 
to the characteristic projected spatial separation between 
<s>- and <L>-scale features, <|A/(r)| 2 > should exhibit a 
minimum, which is indeed the case for scales around 
<d> « 20", i.e., 0.4 pc (see Fig. 11 and Marshak et al. 1997 
for a similar behavior in another physical context). Typical 



Fig. 9. — Wavelet coefficients of Ml-67 at the scale 0'.'3. The orientation 
of the image is the same as of Fig. 1. 


voids of width about 20" are compatible with the appear- 
ance of Ml-67 in Figure 1. 

Since the structure function of order 1 scales as 
< | A f(r) | ) ~ r 0 76 (H x 0.76) at the smallest scales, we infer 
the fractal dimension of Ml-67 to be about 
D = 3 — H x 2.2-2 .3. Thus, Ml-67 appears neither as a 



log 10 (scale) (arcsec) 


Fig. 11. — Structure function analysis of Ml-67. The quantities 
< | A/(r) | p > are plotted down to 3 pixels of WFPC2 (0.3 ”), for p = 1 and 2. 
Three characteristic scales are also indicated (see text). One arcsecond 
corresponds to about 2.2 x 10“ 2 pc assuming a distance of 4.5 kpc for 
Ml-67. 
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smooth, nor space-filling structure. For many different 
nebular objects, we can find estimates of their related fractal 
dimension in the literature. Generally these fractal dimen- 
sions were derived through box-counting techniques on the 
perimeter of the nebula, hence leading to 1 < D p < 2, the 
majority of the values spanning the range 1.2— 1.4 (Bazell & 
Desert 1988; Blacher & Perdang 1990; Falgarone, Phillips, 
& Walker 1991, and references therein). Experiments have 
shown that in the case of isotropic turbulence/fractality, the 
fractal dimension of the perimeter D p and the fractal dimen- 
sion D of the two-dimensional projected nebula satisfy the 
relation: D p = D — 1 (Sreenivasan & Meneveau 1986; 
Beech 1992). Therefore, assuming that the fractal dimension 
D p of Ml-67 does not vary with the angle of projection (i.e., 
assuming isotropy of the turbulent status of the nebula), we 
infer D p (Ml-67) a 1.2-1. 3, which is in good agreement with 
the measurements determined in the above studies and sug- 
gests that this parameter takes on universal values. 

In Figure 12 we show the hierarchy of exponents £(p) as a 
function of p. These exponents have been obtained through 
fitting of the statistical moments over 2.2 decades in scale, 
including the scale break. Note that C(p)/p is no longer con- 
stant when p exceeds 2-2.5. We interpret this fact as evi- 
dence for possible intermittency in the data, i.e., sporadic 
and intense “ bursts ” of high-frequency activity. 

4.2. CFHT Fabry-Perot Data: Structure Function Analysis 
4.2.1. General Results on the Velocity Field 

Figure 13 shows the channels covering the spread in 
radial velocity for the nebula Ml-67. From these data 
Ml-67 appears as a strongly distorted, thick spherical shell 
seen almost exactly along its direction of rapid spatial 
motion in the ISM (Moffat, Lamontagne, & Seggewiss 
1982; Sirianni et al. 1998). The radial velocity of the center 
of expansion is «137 km s -1 (Sirianni et al. 1998). Note 
that towards the star, matter is seen at many different 
velocities. This is due to the broad Ha emission line origin- 
ating in WR 124 (wind terminal velocity: «710 km s -1 ; 
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Fig. 12. — Structure function analysis of Ml-67. The corresponding ((p) 
function demonstrates the multiaffinity of Ml-67. 
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Hamann et al. 1993; Crowther et al. 1995a). Figure 13 also 
shows that the velocity dispersion inside Ml-67 is quite 
large: note that a simple thick shell by itself would not 
explain the multiple radial velocities along the line of sight. 
This velocity dispersion leads one to consider Ml-67 as a 
thick accelerating shell. On the whole, instead of appearing 
as a nice hollow-type shell projected on the sky, we prob- 
ably see the cap of the bow shock nearly straight on from 
behind, the far side being greatly intensity-enhanced (by a 
factor of 8 or more) compared to the near side, probably as 
a result of ramming with the ISM. This was already claimed 
by Solf & Carsenty (1982). 

Despite the presence of a “ protective ” wind-blown cavity 
of radius several 10 pc, expected to have been produced by 
the original O-star progenitor, the current nebula of radius 
~ 1 pc will nevertheless be able to ram the ISM directly. 
This is because of the rapid motion of the central star 
( ~ 200 km s -1 x2 x 10” 4 pc yr _1 ) combined with the 
braking of the cavity expansion in the ISM, allowing the 
nebula to reach the edge of the cavity and enter in contact 
with the ISM. 

The processing of all the interferograms enabled the mea- 
surement of the velocity centroids covering most of the 
surface of the nebula. Because of the low signal-to-noise 
ratio of the approaching, near side of Ml-67, this has been 
done only for the far, “red” side of the nebula, i.e., for 
velocities greater than 137 km s _1 . Figure 14 shows the 
velocity map of Ml-67 for the red component. From this 
figure, Ml-67 exhibits a distorted velocity field reminiscent 
of a thick expanding shell as suggested by Figure 13. In 
addition, the velocity field does not show any clear bipolar 
outflow, contrary to the previous claim by Sirianni et al. 
(1998). We argue that the Sirianni et al. conclusions were 
based on an insufficient amount of data due to poor spatial 
sampling of the velocity field. Note that sudden, huge 
changes in the transparency conditions of the sky occurred 
during the scanning at the end of the night. This explains 
the absence of measurements in the extreme south part of 
Ml-67. 

Figure 15 shows the distribution of the 103,287 measured 
velocity points of the red component. The mean LSR veloc- 
ity is 176.2 km s _1 with a standard deviation of 15.3 km 
s _ 1 . Given the systemic heliocentric radial motion of Ml-67 
(a 137 km s _1 ) and assuming a symmetric expansion of the 
nebula, with a velocity of expansion of about 46 km s _1 
(Sirianni et al. 1998), the maximum expected radial velocity 
should be «183 km s _1 . However, we note numerous 
points above this value (up to about 225 km s _1 ) showing 
the high velocity dispersion within Ml-67. Note that Figure 
15 shows velocities well below 137 km s _1 , down to about 
120 km s _1 , i.e., velocities apparently related to the blue 
component. This is mainly due to velocity perturbations of 
the shell near its periphery, and the resulting difficulty in 
splitting the two components. Moreover, note the presence 
of a small bump near 120-140 km s _1 . Since this value is 
very close to the radial velocity of the center of expansion 
(137 km s _1 ), here we are mainly concerned with regions of 
Ml-67 located near the periphery of the nebula. Such a 
bump is evidence for velocity perturbations likely due to the 
interaction with the interstellar medium. On the whole, the 
irregular nature of the velocity field is possibly related to 
either large variations in the density distribution of the 
ambient ISM, or large variations in the central star mass- 
loss history. 
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Fig. 13. — Maps of the CFHT/SIS-Ha intensity in Ml-67 for the heliocentric radial velocities (in km s -1 ) indicated in the upper right-hand corners. 
Everywhere in the nebula except near the edges, a splitting of the Ha line into at least two components is detected; the high-velocity component is generally 
brighter (by a factor 8 or more) than the low-velocity component. The central star shows up clearly at the same position in each box. 


4.2.2. Structure Function Analysis of the Velocity Field 
Because any valid statistical analysis of velocity fluctua- 
tions has to be done on a (at least local) stationary velocity 
field, one has to subtract global, systematic, nonstationary 
trends from it. Rather than fitting a two-dimensional poly- 
nomial to the velocity data in order to remove the expan- 
sion pattern of the distorted shell (which would generally 
lead to additional new systematic trends), we have con- 


volved the centroid velocity map with Zurflueh filters 
(Zurflueh 1967; Miesch & Bally 1994) of different frequency 
response widths, and removed the smoothed map from the 
original data. The Zurflueh filter is a symmetric, low-pass 
filter which nicely approximates a step function in the fre- 
quency domain and moderates the introduction of artificial 
high-frequency patterns in the convolved data. This system- 
atic trend correction gave us the residual fluctuating veloc- 
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Fig. 14. — Smoothed velocity map of Ml-67 for the heliocentric radial velocities of the “ red ” component (in km s '). North is up, and east is to the left. 
Note the missing chunk of emission in the bottom part due to sudden arrival of clouds near the end of the FP scan. 


ity component on which we have performed the structure 
function analysis. Figure 16 shows the computed second- 
order discrete structure function of the velocity field of 
Ml-67: 

rM - ax + rij J 

xw W) — ' <4) 

where X and r are two-dimensional vectors on the plane of 
the sky, the summation being done on all the pairs separat- 
ed by r, N{r) (Scalo 1984; Miesch & Bally 1994). In Figure 
16 the second-order velocity structure function shows no 



LSR radial velocity (km s -1 ) 

Fig. 15. — Distribution of the 103,287 measured velocity centroids of 
the red component. 


clear, inertial range, as well at small as at larger spatial 
separations. The absence of a well-defined inertial range 
appears very likely because it is independent of the width of 
the applied Zurflueh filter. However, with a distance of 4.5 
kpc to Ml-67 and a pixel size of 0"3, we detect a possible, 
marginally inertial regime occurring in the short range 
16-50 pixels, i.e., 0.1-0.3 pc. A linear least-squares’ fit of the 
slope in the range where there is a possible correlation of 
the results has been done using a power law and gives a 
slope of £(2) » 0.3 (i.e., the related distribution of turbulent 
kinetic energy with spatial scale would be described by an 
energy spectrum, E(k), scaling as k -13 with the wave 
number k). Such a value does not meet the Kolmogorov 
type description for which one would obtain {(2) = 2/3 (i.e., 
E(k) ~ k~ 5/3 ). On the other hand, the 16-50 pixel possible 
inertial range is so short that it is unlikely significant or real. 



LAG (pixels) 

Fig. 16. — Second-order velocity structure function of Ml-67 where 
LAG is the separation between velocity pairs. One pixel corresponds to 
ss6.5 x 10“ 3 pc assuming a distance of 4.5 kpc for Ml-67. 
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At smaller scales, down to the achieved spatial resolution, 
the 2-16 pixel range (i.e., 0.01-0.1 pc) seems to be related 
with an almost stationary, scale-independent field of 
increments. Indeed, in this range S(r) is found to scale 
trivially: ((2) « 0. Note we are not able to confirm the iner- 
tial range detected from the density field (see § 4.1.3) because 
this inertial range ends at about 0'.'3, which is well below the 
resolution of the CFHT data. Similar results are obtained 
when estimating the second-order structure function of the 
centroid velocity field using (linear-trends free) wavelet coef- 
ficients following the method described in § 4.1.1. 

4 . 2 . 3 . The Distribution of Velocity Increments in Ml-67 — 
Structure Function Analysis of the Regularized Velocity Map 

Because we depend on the spatial coordinates to perform 
all averaging operations when analyzing the velocity data, 
our first task is to determine which aspects of the signal are 
more likely to be stationary. The velocity field of Ml-67 is 
nonstationary but with presumably locally stationary 
increments, so that in the previous section we turned to 
structure function analysis. Although possibly related to 
effects of finite spatial resolution (see Marshak et al. 1994, 
and references therein), the small, almost null (at the smal- 


lest scales), value of ((2) derived from the velocity field sug- 
gests that the field of increments Au(A ; r) = v(X) — v(X + r ) 
in equation (4) is apparently essentially scale-independent. 

One-point probability density functions (PDF) make no 
use of the ordering of the data, leaving wide open the 
possibility of correlations between data-values at different 
points or structures in the data. Flowever, the apparent 
scale independence of the velocity increments relaxes this 
drawback of one-point PDF statistics. One traditional 
approach is one-point PDF histograms, typically searching 
for gaussianity. 

Figure 17 shows the empirical PDFs of velocity 
increments as traced by the velocity field wavelet coeffi- 
cients for four different spatial separations. Each PDF is 
compared with a Levy stable distribution (found using the 
computer program STABLE 13 ; Nolan 1997) and the best- 
fit Gaussian distribution. Recall that Levy stable distribu- 
tions arise from the generalization of the central limit 
theorem to a wider class of distributions. In the Gaussian 
case, the renormalized sum of independent and identically 

13 The computer program STABLE can be downloaded from http:// 
academic2.american.edu/~jpnolan/stable/stable.html. 




Fig. 17. — Empirical PDFs of velocity increments as traced by the velocity field wavelet coefficients for four different spatial scale separations (squares). 
Each PDF is compared with a Levy stable distribution (solid curve) and the best-fit Gaussian distribution (dotted curve). The histograms were constructed 
from about 27,300 individual measurements at each scale. 
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distributed random variables with finite variance has the 
same probability distribution as the individual terms in the 
sum. In 1937 Levy showed that the Gaussian distribution is 
but a special case in a wider family of distributions once the 
constraint of finite variance is relaxed (Feller 1971). The 
resulting distributions are referred to as Levy stable dis- 
tributions, or simply Levy distributions, and are character- 
ized by slowly decaying tails and infinite moments. Levy 
distributions, L y , require four parameters to describe: an 
index of stability a in the interval (0, 2], a skewness param- 
eter — 1 < ft < 1, a scale parameter a > 0, and a location 
parameter /z. The index a determines the rate at which the 
tails of the distribution decrease [L a (x) oc x _1_a ]; when 
a = 2, a Gaussian distribution results. When a < 2 the 
variance is infinite. In general, the pth moment of a Levy 
stable random variable is finite if and only if p < a. The 
parameters a and /z determine the width and the shift of the 
peak of the distribution. For a positive (respectively, 
negative) p the distribution is skewed to the right 
(respectively, the left). For ft = 0, the distribution is sym- 
metric about /z. In addition, as a approaches 2, [i loses its 
effect and the distribution approaches a Gaussian regard- 
less of fi (Nolan 1997). 

Strictly speaking, velocity increments cannot have a Levy 
stable distribution because this would imply nonzero prob- 
ability density in regions outside the range for which veloc- 
ity increments are physically meaningful. However, these 
distributions can yield useful approximations over a 
truncated/finite range as demonstrated in other physical 
contexts (Painter, Berersford, & Paterson 1995). In addi- 
tion, when we use an infinite variance distribution in order 
to fit the distribution of velocity increments, all that we are 
saying is that the observed slowly decaying tails result in 
behavior similar to that of the infinite variance model. On 
the whole, we stress the fact that variance is but one 
measure of spread for a distribution, and is not appropriate 
for all problems (Painter et al. 1995). 

As the scale increases, Figure 17 shows us that the dis- 
tribution approaches a Gaussian distribution, as expected 
for a truncated Levy distribution (Nakao 2000). At small 
scales, the distribution of the velocity increments is charac- 
terized by a nonnegligible number of extreme events 
making the PDF of the velocity increments not narrow 
enough to enable a simple dimensional argument to relate 
quantitatively all moments : < | Azz(r) | p > no longer approx- 
imates < | Azz(r) | y and precludes the determination of any 
actual (multi)fractal structure. At large scales the gauss- 
ianity of the PDFs releases this drawback and nontrivial 
scaling laws may be detected, although likely biased, as seen 
in Figure 16. Note that the Levy fits in Figure 17 are nearly 
symmetric (j? « 0) at small scales. However, as the scale 
increases, ji increases mainly because the wavelet coeffi- 
cients take into account more and more nonlinear trends in 
the signal. 

Because of the high spatial resolution of our Fabry-Perot 
observations and the good seeing conditions, one expects 
each pixel to look at a fairly independent piece of the sky. 
However, very low flux regions could see their spectrum 
polluted by the point-spread function wings of a very bright 
neighbor. Such combination is indeed found in our data 
and leads to the dominant, larger velocity increments, in 
line with the <s>-scale intense features detected in the 
density field wavelet analysis (see § 4.1.3). As a consequence, 
the velocity centroid measurements are sparsely corrupted 


resulting in fat-tail PDFs of velocity increments. Note that, 
although there is a regular grid of pixels, the emission is 
distributed irregularly according to the density distribution. 
Should the velocity signal be coming primarily from the 
brightest regions then the scaling signal could be addi- 
tionally confused by the irregular grid signal (see Rauzy, 
Lachieze-Rey, & Henriksenet al. 1995). However, our FP, 
high throughput velocity map is not found to suffer from a 
significant interference signal from a nonuniform grid. 

We therefore applied a 1 pixel radius ring median filter to 
the velocity map. The ring median filter is defined as a 
median filter that assigns weights only to selected pixels in 
an annulus. Its advantage is that it has a sharply defined 
scale length; that is, all objects with a scale size less than the 
radius of the ring are filtered and replaced by the local 
background level. It provides a simple method to remove 
noiselike deviations (independent of morphology) from a 
digital image, leaving behind the large-scale structures and 
overall gradients (Seeker 1995). We expect each velocity 
point of the regularized velocity map to be an accurate 
measure of the underlying gas kinematics, averaged on its 
pixel seeing disk. 

Figure 18 shows the computed structure functions of the 
velocity field of Ml-67 after small-scale regularization using 
a 1 pixel radius ring median filter. The scaling range (a 1"4, 
10") where the exponents are defined is clearly marked. The 
break that defines the upper limit of «10" is not robust 
with respect to adding to the average more data points, and 
simply reflects a violation of the ergodicity assumption. 
Below the 174 scale, the slope is slightly steeper, as a conse- 
quence of the ring median filtering, which tends to smooth 
the signal at the smallest scales. Indeed, using a 2 pixel ring 
annulus simply results in moving the lower limit for scaling 
towards larger scales. Therefore, the cleaned velocity map is 



Fig. 18. — Velocity structure function analysis of Ml-67 after small- 
scale regularization using a 1 pixel radius ring median filter. The quantities 
< | At >(r) | p > are plotted down to 3 pixels of CFHT/SIS (0.9”), for p = 0.5- 
2.5, in steps of 0.5. One arcsecond corresponds to about 2.2 x 10 -2 pc 
assuming a distance of 4.5 kpc for Ml-67. 
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found to scale over at least 1 order of magnitude in scale, 
with ((1) « 0.48-0.49 and £(2) « 0.90-0.91. The latter value 
is close to 1.0 (E(k) ~ k~ 19 ), suggesting the dynamics is 
essentially shock dominated (for which E(k) ~ k~ 2 ). Figure 
19 shows the function £ as a function of p, up to 8.0, in steps 
of 0.1. Clearly £(p)/p is nonlinear and suggests the multi- 
affinity of the velocity field. 

Multiaffine modeling could be done in the context of 
“ universal multifractals ” (Schertzer & Lovejoy 1987) where 
multifractals are the generic result of multiplicative cas- 
cades. A continuous-scale limit of such processes leads to 
the family of log-infinitely divisible distributions that have a 
Gaussian or Levy stable generator (Schertzer & Lovejoy 
1987, 1997; Lovejoy & Schertzer 1990). For universal multi- 
fractals we have 

£(p)=P«l)--^-(P‘-p)(«/l), (5) 

a — 1 

C(p) = K(l) - Ci p In (p) (a = 1) , (6) 

where C l < 2 is an intermittency parameter and 0 < a < 2 
is the Levy tail index. After the existence of a scaling regime 
is established, one can test for universal multifractal behav- 
ior. An efficient technique for that purpose is the “ double 
trace moment” (DTM) method (Lavallee 1991). A new 
function K(q, q) is first defined as K(qq) — qK(q), with 
K(q) = qC( 1) — £(< 2 ). In the case of universal multifractals, 
the ^-dependence in K(q, q) factorizes as K(q, q) = q'K(q). 
The DTM method allows the computation of K(q, q), and 
assuming universality, the Levy index a can be estimated by 
fixing q and varying q. With a good estimate of a in hand, an 
ordinary least-squares estimate of C 1 is easy to obtain. In 
our case, a standard two-parameter nonlinear regression 
does very fine over almost 2 orders of magnitude in q, with 



p 


Fig. 19. — Velocity structure function analysis of Ml-67 after small- 
scale regularization using a 1 pixel radius ring median filter. The corre- 
sponding {(p) function (points ) demonstrates the multiaffinity of the 
velocity field of Ml-67. A universal multifractal fit is also shown as a solid 
curve (see text). 


lO 



Fig. 20. — Plots of log 10 \K(q, q)\ (using q = 0.5, 1.5, 2, 2.5, 3, and 5.5) 
obtained with a DTM analysis of the velocity field of Ml-67 after small- 
scale regularization. From the slopes and the log 10 q = 0 intercepts of the 
linear region we obtain a = 1.90-1.92, and C, = 0.04 + 0.01. There is a 
first-order multifractal phase transition near max (qq, q) = 7, which is a 
consequence of the finiteness of the analyzed sample. 

a « 1.91 and C, « 0.04 + 0.01 (see Fig. 20). The related £(p) 
universal multifractal fit is shown in Figure 19 ( solid curve). 
The agreement between the empirical £(p) and the multi- 
fractal fit is very good for p below «7. Note that the finite 
size of the sample implies that sufficiently high-order 
moments are dominated by the largest values in the field 
and therefore underestimate the true ensemble moments. 
Beyond max (qq, q) = (2/C, l 1 '* « 7.7 equation (5) is no 
longer true, and £(p) becomes linear (see Fig. 19). This is the 
so-called first-order multifractal transition (Schertzer & 
Lovejoy 1992), which is also seen in Figure 20 where the 
empirical K(q, q fs are no longer fitted by the linear regres- 
sions. 

5. DISCUSSION AND CONCLUSION 

Using the unique, high-resolution, wide-field imaging 
capabilities of HST, we have tested quantitatively, via 
scaling laws, for compressible turbulence in the distribution 
of density clumps in the ejected nebula Ml-67. These data 
have been straddled by complementary CFHT Fabry-Perot 
Ha data in order to perform the same test on the velocity 
field of Ml-67. What is generally understood as turbulence 
is the difference resulting from the subtraction of the instru- 
mental and thermal broadening from the observed emission 
line widths. In other words, turbulence is a quantity related 
to macroscopic, chaotic motions of the gas and their related 
density fluctuations. 

At small scales, the density field of Ml-67 appears 
remarkably structured in essentially chaotically oriented fil- 
aments everywhere in the nebula. However, many filaments 
(especially at the smallest scales) seem to be slightly radially 
oriented. This is not surprising, since the force that most 
likely drives the nebula and the ensuing turbulence is 
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directed outwards in all directions from the central star. The 
presence of filaments agrees with earlier numerical simula- 
tions of turbulent flows (e.g., Sreenivasan & Antonia 1997; 
Vazquez-Semadeni 1999, and references therein) or obser- 
vations (e.g., Miville-Deschenes et al. 1999, and references 
therein), which have revealed how turbulence is character- 
ized by filamentary vortical structures and surrounding dis- 
sipative sheets. It is worth noting that the empirical PDFs 
in Figure 17 are statistically good estimates of the true 
PDFs in the central parts of the distributions. Clearly, these 
PDFs are not Gaussian at the smallest scales; rather, they 
are best described by Levy stable distributions that are 
characterized by extended tails. In addition, as the spatial 
scale increases, the Levy index a is found to increase, the 
Gaussian case (a = 2) being attained for scales above about 
4". This behavior also suggests the presence of vorticity in 
the flow and is likely related to the phenomenon of inter- 
mittency (Lis et al. 1998, and references therein). 

The velocity field of Ml -67 (after correction for bad 
velocity points related to the largest velocity increments) 
shows a clear inertial range reminiscent of turbulence. 
Because the ISM and the ejected nebula are compressible 
and inhomogeneous, the gas flows being in addition highly 
supersonic, we expect shock wave propagations within 
Ml -67. External forces due to stellar winds (or gravitation 
in another astrophysical context, such as giant molecular 
clouds: Gill & Henriksen 1990) play an important role in 
feeding the ISM with highly pressurized motions, moving 
the gas at velocities larger than the sound speed, and ulti- 
mately leading to supersonic turbulence and its related 
small-scale compressions. Interaction between shock waves 
and turbulent fluctuations are expected to distribute the 
energy dissipation on a cascade characterized by a steeper 
law C(l) = 0.5 (Fleck 1996) than predicted by Kolmogorov’s 
law K(l) = |]. Our value of £(1) = 0.48-0.49 strongly sug- 
gests that the turbulence traced by the velocity field is com- 
pressible. In addition, since our observations come from a 
projection on the plane of the sky, we therefore expect a 
smearing of the information along each line of sight. This 
may lead to additional perturbations in the flows, with 
apparently fluctuating power laws as a function of projected 
spatial separation (O’Dell & Castaneda 1987, and refer- 
ences therein) precluding any inertial regime detection. 
Figure 18 tells us that this effect is practically negligible, and 
additionally suggests that the Ha layer is thin compared to 
the projected spatial separations. On the other hand, the 
value £(2) = 0.90 is in very good agreement with the empiri- 
cal Larson-type law, which predicts ((2) = 0.8 ± 0.1 
(Larson 1981; Miesch & Bally 1994) and is associated to 
fields that are essentially shock-dominated. Note that 
C(2) < 2 x C(l) results from the slight intermittency of the 
velocity field, which makes it multifractal rather than 
monofractal. With C 1 close to 0, the process is found to be 
nearly homogeneous (the fractal dimension of the set con- 
tributing to the mean velocity field is 2 — C 1 x 1.96). On 
the other hand, the degree of multifractality a « 1.90-1.92 is 
close to 2. The latter result suggests that we are far from the 
so-called monofractal /? -model (a = 0); rather, the turbu- 
lence is best described by the log-normal model (a = 2) for 
which high values of the field dominate more than for 
smaller values of a. The estimation of a and C 1 for other H ii 
regions appears now essential in order to test whether the 


multifractal parameters found in the case of Ml -67 take on 
universal values. For comparison with other turbulence 
studies, one may derive from C 1 the standard intermittency 
parameter n which is the autocorrelation exponent of the 
dissipation field: n = K( 2,1). For the /1-model (a = 0), n = 
C u whereas for the lognormal model (a = 2), n = 2 x C 1 . 
In our case, with a « 1.91, we get n « 1.90 x C, « 0.076. 

Recall that turbulent velocity measurements in wind 
tunnel experiments lead to a Kolmogorov-type (i.e., incom- 
pressible; £(2) = 2/3) turbulence with a x 1.3 ± 0.1, C x x 
0.25 ± 0.05 and n x 0.35 ± 0.1 (Schmitt et al. 1992). The 
higher degree of multifractality and the smaller inter- 
mittency found for the velocity field of Ml -67 is likely a 
consequence of the compressible nature of the tracing gas. 
This will be discussed in a forthcoming paper. 

It is worthy of note that large variations may exist in the 
history of WR 124’s mass-loss and the distribution of the 
ambient ISM. This precludes (1) any unique interpretation 
of the absence of a well-defined inertial range in the density 
field and (2) any unambiguous derivation of the stellar 
mass-loss history because density fluctuations in the 
absence of self-gravitational effects are believed to be tran- 
sient only in the context of supersonic turbulence (Vazquez- 
Semadeni 1994). However, note that the marginally inertial 
range found from the density field is likely a consequence of 
the rapid regeneration of density structures triggered by the 
turbulent nature of the velocity field. Finally, given the 
extreme perturbations of the velocity and density fields in 
Ml-67, it is virtually impossible to measure any systematic 
impact of the WR (or LBV) wind on the nebular structure. 
For whatever reason, any systematic effects seem to have 
been randomized such that they have become fully masked. 
On the other hand, the origin of the <s)-scale features found 
from the wavelet analysis of the HST Ha image still has to 
be investigated in detail. These features could be related to 
either clumps of stellar origin, or smoothed out versions of 
the bright “ bullets ” reported in Grosdidier et al. (1998). On 
the other hand, the irregular nature of the velocity field is 
likely due to either large variations in the density distribu- 
tion of the ambient ISM, or large variations in the central 
star mass-loss history. 

On the whole, Ml-67 exhibits the signatures of compress- 
ible, intermittent turbulence, the driver being ultimately the 
stellar wind originating in WR 124, as opposed to gravity/ 
magnetic fields in the ISM. In addition, from the structure 
function of order 1 of the density field (which scales as 
< | A/(r) | ) ~ r 0 - 76 at the smallest scales), we infer the fractal 
dimension of Ml-67 to be D x 2.2-23. This result is very 
close to the values often quoted for molecular clouds and 
indeed is essentially the result found in the wavelet analysis 
of Gill & Henriksen (1990). Therefore, although different 
drivers may be related to the generation of supersonic, com- 
pressible turbulence, it appears that the resulting phenom- 
enology essentially remains the same. 
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